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We investigate the dynamics of two models of biological networks with purely suppressive in¬ 
teractions between the units; species interacting via niche competition and neurons via inhibitory 
synaptic coupling. In both of these cases, power-law scaling of the density of states with probability 
arises without any fine-tuning of the model parameters. These results argue against the increasingly 
popular notion that non-equilibrium living systems operate at special critical points, driven by there 
by evolution so as to enable adaptive processing of input data. 

PACS numbers: 


Are biological systems adjusted by evolution to oper¬ 
ate at a special critical point? [T] This idea has surfaced 
in variety of different contexts over the past years, in¬ 
cluding network architectures scaling in neural 

processing m n E], and self-organization of biological 
flocks laiT]. Underlying this concept is the appealing but 
vague notion that being critical is somehow conducive to 
flexible adaptation in the face of large variations in stim¬ 
uli to be acted upon by the system in question. 

One approach to demonstrating this special critical na¬ 
ture relies upon a maximum entropy approach to fit cor¬ 
relations of the fundamental events miHi, say spikes in 
the neural context. This method constructs an effective 
Hamiltonian for the system that reproduces the equal¬ 
time correlations and expectation values. This Hamilto¬ 
nian then allows for the evaluation of the specific heat as 
a function of temperature T, where by definition T = 1 
is used to fit the original data. A specific heat peak at 
T = 1 which grows as a function of system size is then in¬ 
terpreted as evidence of criticality of the actual biological 
process. 

Here, we explore the extent to which this picture is 
valid, using two different classes of biological models aris¬ 
ing respectively in the ecological and neural contexts. 
These models have in common that the nonlinear inter¬ 
action between the fundamental entities are purely in¬ 
hibitory. Species inhibit each other’s survival by com¬ 
peting for the same niche [9] and also our neurons have 
purely inhibitory couplings as they attempt to classify 
incoming activation signals. In both of these cases, non¬ 
trivial power-law scaling of the integrated density of 
states Af{p) with the probability p emerges without hav¬ 
ing to carefully adjust parameters or consider sophisti¬ 
cated feedback scenarios. In the language of Ref. [8], the 
entropy S = logM{p) is linear in the energy. Our sys¬ 
tems are inherently non-equilibrium, with a power-law 
distribution of net fluxes (see later) and hence cannot 
be described by any Hamiltonian. Forcing an equilib¬ 
rium Hamiltonian to reproduce the expectation values 
and equal-time correlations leads to a theory that looks 


critical in the aforementioned sense, but this feature is 
not a useful characterization of the actual operating state 
of these processes. 

Generalized Competitive Lotka-Volterra Sys¬ 
tem We start by reviewing a recent study m of a gen¬ 
eralized Lotka-Volterra model m to study competing 
ecological species mm- This model describes a semi- 
isolated island with infrequent immigration of individuals 
from an external “mainland” on which reside Q species. 
The individuals on the island undergo a stochastic birth- 
death process, where the competition-induced death rate, 
di, for individuals of a given species i is determined by 
the abundances, n/., /c = 1,... Q of all the various species 
(including its own) on the island, 

= ^ + ’ ( 1 ) 

where K represents the carrying capacity and C measures 
the strength of competition. The interaction strengths 
Cij are chosen randomly from a Gamma distribution with 
mean unity and standard deviation a. For simplicity, 
the birth rates of all the different species are chosen to 
be identical, a value which is set equal to unity to de¬ 
fine the time scale. The immigration rates of the various 
species are also fixed to be identical, and denoted by 
A. This model was studied intensively m, and found 
to exhibit four different regimes, or “phases”, of behav¬ 
ior as the competition strength C is varied. For small 
C, the competition effects are minor and typically all 
Q species are present on the island in steady-state (the 
full coexistence phase). At some critical C, the weakest 
species basically goes extinct, except for occasional ulti¬ 
mately unsuccessful attempts to recolonize the island via 
immigration. Increasing C further causes more species 
to undergo extinction (the partial coexistence phase). 
At a second critical C, the dynamics changes dramat¬ 
ically, and there is no longer any set of species which 
is permanently represented. Rather, the system under¬ 
goes regime shifts where, due to the combined effects of 
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FIG. 1: Snapshots of the abundance for all species in the 
stochastic model for various levels of competition, C. Red is 
high abundance, dark blue low, as shown in the color bar at 
the right. Q = 20, iF = 100, a = 0.5, A = 0.01. The time is 
measured in units of 0.08/A. The first panel exhibits full co¬ 
existence, the second and third partial coexistence, the fourth 
and fifth disorder dynamics, and the last glass-like dynamics. 
(From Ref. [TO]) 


extinction and immigration, the set of dominant species 
changes completely from time to time (the disordered 
phase). Lastly, at C’s close to unity, the system is marked 
by long periods where some small set of species coexists, 
with infrequent regime shifts between a few such sets (the 
glass-like phase). This behavior is illustrated in Fig. 
where for various C’s, snapshots of the abundance of all 
Q = 20 species through time is presented. 

In order to meaningfully measure the entropy, we need 
to reduce the dimensionality of the state space. Essen¬ 
tially, we wish to classify a state by which species are 
stably present. To do this it is necessary to eliminate the 
very low abundance species, with only 1 or 2 representa¬ 
tives, since their presence is in general very transitory. A 
convenient way to accomplish this without introducing an 
arbitrary threshold is via the instantaneous “quenching” 
procedure described in Ref. [10]. Each snapshot is taken 
as the initial condition of a deterministic run with A = 0, 
with equations of motion hi = {l — di)ni^ which is run un¬ 
til it reaches a steady-state. In this way, low abundance 
species with negative net growth rates are eliminated. We 
then calculate the number of appearances of every state 
in a long run, which when normalized by the number of 
snapshots gives the probability p of that state. Erom this, 
we generate the complementary cumulative distribution 
function, A/’(p), the number of states whose probability 
is greater than a given p. The results of a typical run for 
Q = 40 is given in Eig. We see that Af{p) behaves as 
a power-law for all but the largest p’s, over a span of 4.5 
decades. This translates directly to an asymptotically 
linear S{E) relation as E = — logp ^ +oo. As discussed 
in detail in Ref. [5] , this means that the attempt to find 
a typical value of E at the fixed temperature T = 1 via 


FIG. 2: The complementary cumulative distribution function 
JV{p); i.e., the number of states with probability greater than 
p, for runs of 10® or 10^ consecutive snapshots, presented 
in log-log scale. Q = 40, A = 200, C = 0.35, a = l/\/8, 
A = 0.01. The time between snapshots is 0.08/A. Also shown 
is a power-law fit to the low-p data. 
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FIG. 3: The complementary cumulative distribution func¬ 
tion J\f{p)] i.e., the number of states with probability greater 
than p, for C = 0.25,0.45,0.55 with Q = 60 in (a) and 
Q — 20,40,60 with C — 0.55 in (b). The other parameters 
are A = 0.01, K — 50Q, a — .5y/20/Q. The time between 
snapshots is 0.08/A. 

the thermodynamic relation = 1 cannot succeed, and 
the system exhibits a rather unusual type of critical be¬ 
havior. 

One question that immediately arises is how general 
this phenomenon is, or whether it holds only for a spe¬ 
cific set of parameters or range of parameters. To inves¬ 
tigate this, we varied C in the range from 0.35 to 0.55, 
measuring M{p) for each. We also varied Q, the number 
of interacting species. As Q is changed, we scale the vari¬ 
ance parameter a of the inhibition matrix as Ijy/Q [9| 
and K linearly with Q. The results are presented in Fig. 

Of interest is the fact that the linearity improves as 
Q increases and that the data approaches a pure Zipf 
law [13 HU with scaling exponent unity. The fact that 
“criticality” extends for a finite parameter range in C 
is presumably due to the inherently non-equilibrium na¬ 
ture of the problem. One way to measure deviations from 
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FIG. 4: Distribution of net fluxes (in units of the immigration 
rate A) for the case Q = 20, C = 0.45, a = 0.5, A = 0.01. 

equilibrium is via the violation of detailed balance in the 
flux between states. In Fig. we present a plot of the net 
steady-state flux between different states in the model, 
showing that these again have a power-law distribution. 
Of course these would all be precisely zero in any model 
where the steady-state obeyed detailed balance. 

Neural Net Model The main nonlinear ingredient of 
our ecology model is the dynamic competition between 
different species. We are thus motivated to turn to a 
model with a similar interaction in a different context, 
a set of neurons which are mutually inhibitory and are 
driven by external sources [15] . The inhibition matrix Jij 
is constructed from random elements of average size unity 
(distributed uniformly between 0 and 2), with density 
p, with all other elements set to zero. As in standard 
integrate and fire networks, the voltage level on a given 
neuron decays exponentially with time constant r. The 
external input, V/ is uniformly distributed at each time 
slice dt between (it(l± v^cr/). The threshold for a neuron 
firing is given initially by 

Vc=t(i- ( 2 ) 

SO that in the absence of noise, every neuron would fire 
with period Tp. Thus, the update rule for Vi is 

Vi{t + dt) = + v/ (3) 

At this point, suprathreshold neurons fire, and are reset 
to 0. The firing neurons then reduce the voltage on the 
neurons they are connected to according to 

Vi = Vi-Cj2 Jij (4) 

j efire 

A record is kept of every neuron that fires during each 
window of duration Tw, which in this case is reasonable 
to take as T^, since this is the minimum interval between 
firings (modulo the input noise). 


A pictorial view of the dynamics is given in Fig. 
Notice that the dynamical state of the system is very 
different than that depicted in Fig Ic, which presents a 
similar trace for the colony model in the range of param¬ 
eters we have been considering. Essentially, the strong 
inhibition and the lack of any long-term memory that 
a neuron has been active serves to create an extremely 
chaotic process; in the ecology model the role of mem¬ 
ory is played by the size of a species population which if 
large prevents the species from rapidly disappearing. To 
mimic this effect, we added a feedback between firing and 
threshold. When a neuron fires, its threshold for future 
firing is lowered, by a factor /t- The threshold then re¬ 
laxes exponentially back to its nominal value, Vc, with a 
relaxation rate of 7 ^. One could imagine accomplishing 
analogous changes by altering synapses, i.e. weakening 
the synaptic inhibition for any neuron that fires. Once 
this change is implemented, the dynamics stabilizes to 
the form seen in Fig. dt- Somewhat remarkably, the 
actual variation in thresholds achieved by our algorithm 
once the system has settled into a statistical steady-state 
is only about 10% (data not shown). Yet, this is enough 
to create robust patterns of firing and non-firing neurons 
which are then acted upon by the fluctuation input sig¬ 
nals, since the period of firing of a neuron is exponentially 
sensitive to its threshold when the threshold is close to its 
critical value of r. This capability of using relatively weak 
induced memory to create stable firing patterns should 
be studied in more detail in future works. 

We then tabulate the statistics of firing patterns. Our 
first set of runs is with = 80 neurons and parameters 
p = 1, T = 20, Tp = 83, C = 1, a/ = 0.02. In Fig. we 
again turn to the complementary cumulative distribution 
function. As in the previous model, we now observe an 
almost perfect Zipf law form [16] , scaling with exponent 
close to -1. This would guarantee that any attempt to 
model the thermodynamics of the system with an equilib¬ 
rium Hamiltonian would “discover” that the system looks 
like it had been tuned to be very close to a critical point. 
Again however, this state exists over a broad range of 
parameters and does not require any fine-tuning. A sim¬ 
ilar tabulation for the non-adaptive model shows that 
the system thermodynamics is completely dominated by 
a huge number of very low occupancy states. Approaches 
which claim that scaling should exist in most models of 
this kind [14] need to be able to explain what fails in our 
non-adaptive neural system. 

Scaling behavior of various kinds seems to be a ubiq¬ 
uitous feature of many non-equilibrium systems [17] . We 
have seen here two examples of such systems. We are 
certainly not claiming that all or even most models one 
could write down will automatically have such a non¬ 
trivial property. As shown in our neural example, we 
need to have the proper mix of memory and transitions; 
too little memory and the system has no enduring high 
probability states and too much memory (for example, 
by lowering the inhibition) locks the system permanently 
into too few states. But unlike the expectation for equi- 
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FIG. 5: Sample firing patterns from our neural net model with 
N = 80 neurons. Left: Pattern with constant threshold Vc = 
19.67. The time window was Tw = Tp = 83. Right: Pattern 
with dynamical thresholds, governed by the parameters fr = 
0.85, 7t = 0.02. Due to the lower threshold of the adapted 
neurons, the time window was taken as Tw = 40. For both 
models, C = 1, a = 0.7, aj = 0.02, r = 20. 


A^=80 



FIG. 6: Gomplementary cumulative distribution function for 
N = 80. Parameters are the same as in the simulation shown 
on the right side of Fig. 


librium distributions, there is a ’’Goldilocks” range, not 
a ’’Goldilocks” point. This of course would make it much 
easier for evolution to drive biological processes to have 
these properties, should they prove advantageous for in¬ 
formation processing in response to stimuli. 
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